module GOLD_domains

!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of GOLD.                                        *
!*                                                                     *
!* GOLD is free software; you can redistribute it and/or modify it and *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* GOLD is distributed in the hope that it will be useful, but WITHOUT *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

#undef  FORBID_NONBLOCKING_UPDATES

use GOLD_error_handler, only : GOLD_error, NOTE, WARNING, FATAL, is_root_pe
use GOLD_file_parser, only : read_param, get_param, log_param, log_version
use GOLD_file_parser, only : param_file_type

use GOLD_coms, only : PE_here, root_PE, num_PEs, GOLD_infra_init, GOLD_infra_end
use GOLD_coms, only : broadcast, sum_across_PEs, min_across_PEs, max_across_PEs
use mpp_domains_mod, only : mpp_define_layout
use mpp_domains_mod, only : GOLD_define_io_domain => mpp_define_io_domain
use mpp_domains_mod, only : GOLD_define_domain => mpp_define_domains
use mpp_domains_mod, only : domain2D, domain1D, mpp_get_data_domain
use mpp_domains_mod, only : mpp_get_compute_domain, mpp_get_global_domain
use mpp_domains_mod, only : global_field_sum => mpp_global_sum
use mpp_domains_mod, only : mpp_update_domains, CYCLIC_GLOBAL_DOMAIN, FOLD_NORTH_EDGE
#ifndef FORBID_NONBLOCKING_UPDATES
use mpp_domains_mod, only : mpp_start_update_domains, mpp_complete_update_domains
#endif
use mpp_parameter_mod, only : AGRID, BGRID_NE, CGRID_NE, SCALAR_PAIR, BITWISE_EXACT_SUM, CORNER
use mpp_parameter_mod, only : To_East => WUPDATE, To_West => EUPDATE
use mpp_parameter_mod, only : To_North => SUPDATE, To_South => NUPDATE

implicit none ; private

#include <GOLD_memory.h>

public :: GOLD_domains_init, GOLD_infra_init, GOLD_infra_end, get_domain_extent
public :: GOLD_define_domain, GOLD_define_io_domain
public :: pass_var, pass_vector, broadcast, PE_here, root_PE, num_PEs
public :: pass_var_start, pass_var_complete
public :: pass_vector_start, pass_vector_complete
public :: global_field_sum, sum_across_PEs, min_across_PEs, max_across_PEs
public :: AGRID, BGRID_NE, CGRID_NE, SCALAR_PAIR, BITWISE_EXACT_SUM, CORNER
public :: To_East, To_West, To_North, To_South, To_All

interface pass_var
  module procedure pass_var_3d, pass_var_2d
end interface pass_var

interface pass_vector
  module procedure pass_vector_3d, pass_vector_2d
end interface pass_vector

interface pass_var_start
  module procedure pass_var_start_3d, pass_var_start_2d
end interface pass_var_start

interface pass_var_complete
  module procedure pass_var_complete_3d, pass_var_complete_2d
end interface pass_var_complete

interface pass_vector_start
  module procedure pass_vector_start_3d, pass_vector_start_2d
end interface pass_vector_start

interface pass_vector_complete
  module procedure pass_vector_complete_3d, pass_vector_complete_2d
end interface pass_vector_complete

type, public :: GOLD_domain_type
  type(domain2D), pointer :: mpp_domain => NULL() ! The domain with halos on
                                        ! this processor, centered at h points.
  integer :: nxtot, nytot               ! total horizontal domain sizes
  integer :: nx_halo, ny_halo           ! X- and Y- halo sizes in memory.
  logical :: symmetric                  ! True if symmetric memory is used with
                                        ! this domain.
  logical :: nonblocking_updates        ! If true, non-blocking halo updates are
                                        ! allowed.  The default is .false. (for now).
  integer :: layout(2), io_layout(2)    ! Saved data for sake of constructing
  integer :: X_FLAGS,Y_FLAGS            ! new domains of different resolution.
  logical :: use_io_layout              ! True if an I/O layout is available.

end type GOLD_domain_type

integer, parameter :: To_All = To_East + To_West + To_North + To_South

contains

subroutine pass_var_3d(array, GOLD_dom, sideflag, complete, position)
  real, intent(inout), dimension(:,:,:) :: array
  type(GOLD_domain_type), intent(inout) :: GOLD_dom
  integer, optional, intent(in) :: sideflag
  logical, optional, intent(in) :: complete
  integer, optional, intent(in) :: position
! Arguments: array - The array which is having its halos points exchanged.
!  (in)      GOLD_dom - The GOLD_domain_type containing the mpp_domain needed to
!                      determine where data should be sent.
!  (in)      sideflag - An optional integer indicating which directions the
!                       data should be sent.  It is TO_ALL or the sum of any of
!                       TO_EAST, TO_WEST, TO_NORTH, and TO_SOUTH.  For example,
!                       TO_EAST sends the data to the processor to the east, so
!                       the halos on the western side are filled.  TO_ALL is
!                       the default if sideflag is omitted.
!  (in)      complete - An optional argument indicating whether the halo updates
!                       should be completed before progress resumes.  Omitting
!                       complete is the same as setting complete to .true.
!  (in)      position - An optional argument indicating the position.  This is
!                       usally CORNER, but is CENTER by default.
  integer :: dirflag
  logical :: block_til_complete

  dirflag = To_All ! 60
  if (PRESENT(sideflag)) then ; if (sideflag > 0) dirflag = sideflag ; endif
  block_til_complete = .true.
  if (present(complete)) block_til_complete = complete

  call mpp_update_domains(array, GOLD_dom%mpp_domain, flags=dirflag, &
                          complete=block_til_complete, position=position)

end subroutine pass_var_3d


subroutine pass_var_2d(array, GOLD_dom, sideflag, complete, position)
  real, intent(inout), dimension(:,:) :: array
  type(GOLD_domain_type), intent(inout) :: GOLD_dom
  integer, optional, intent(in) :: sideflag
  logical, optional, intent(in) :: complete
  integer, optional, intent(in) :: position
! Arguments: array - The array which is having its halos points exchanged.
!  (in)      GOLD_dom - The GOLD_domain_type containing the mpp_domain needed to
!                      determine where data should be sent.
!  (in)      sideflag - An optional integer indicating which directions the
!                       data should be sent.  It is TO_ALL or the sum of any of
!                       TO_EAST, TO_WEST, TO_NORTH, and TO_SOUTH.  For example,
!                       TO_EAST sends the data to the processor to the east, so
!                       the halos on the western side are filled.  TO_ALL is
!                       the default if sideflag is omitted.
!  (in)      complete - An optional argument indicating whether the halo updates
!                       should be completed before progress resumes.  Omitting
!                       complete is the same as setting complete to .true.
!  (in)      position - An optional argument indicating the position.  This is
!                       usally CORNER, but is CENTER by default.

  integer :: dirflag
  logical :: block_til_complete

  dirflag = To_All ! 60
  if (PRESENT(sideflag)) then ; if (sideflag > 0) dirflag = sideflag ; endif
  block_til_complete = .true.
  if (present(complete)) block_til_complete = complete

  call mpp_update_domains(array, GOLD_dom%mpp_domain, flags=dirflag, &
                        complete=block_til_complete, position=position)

end subroutine pass_var_2d

function pass_var_start_2d(array, GOLD_dom, sideflag, position, complete)
  real, dimension(:,:),   intent(inout) :: array
  type(GOLD_domain_type), intent(inout) :: GOLD_dom
  integer,      optional, intent(in)    :: sideflag
  integer,      optional, intent(in)    :: position
  logical,      optional, intent(in)    :: complete
  integer :: pass_var_start_2d
! Arguments: array - The array which is having its halos points exchanged.
!  (in)      GOLD_dom - The GOLD_domain_type containing the mpp_domain needed to
!                      determine where data should be sent.
!  (in)      sideflag - An optional integer indicating which directions the
!                       data should be sent.  It is TO_ALL or the sum of any of
!                       TO_EAST, TO_WEST, TO_NORTH, and TO_SOUTH.  For example,
!                       TO_EAST sends the data to the processor to the east, so
!                       the halos on the western side are filled.  TO_ALL is
!                       the default if sideflag is omitted.
!  (in)      position - An optional argument indicating the position.  This is
!                       may be CORNER, but is CENTER by default.
!  (in)      complete - An optional argument indicating whether the halo updates
!                       should be initiated immediately or wait for second 
!                       pass_..._start call.  Omitting complete is the same as 
!                       setting complete to .true.
!  (return value) - The integer index for this update.
  integer :: dirflag

  dirflag = To_All ! 60
  if (PRESENT(sideflag)) then ; if (sideflag > 0) dirflag = sideflag ; endif

#ifdef FORBID_NONBLOCKING_UPDATES
  call GOLD_error(FATAL, "pass_var_start: This call is not yet available.")
  pass_var_start_2d = -1
#else
  pass_var_start_2d = mpp_start_update_domains(array, GOLD_dom%mpp_domain, &
                          flags=dirflag, position=position)
#endif
end function pass_var_start_2d

function pass_var_start_3d(array, GOLD_dom, sideflag, position, complete)
  real, dimension(:,:,:), intent(inout) :: array
  type(GOLD_domain_type), intent(inout) :: GOLD_dom
  integer,      optional, intent(in)    :: sideflag
  integer,      optional, intent(in)    :: position
  logical,      optional, intent(in)    :: complete
  integer                               :: pass_var_start_3d
! Arguments: array - The array which is having its halos points exchanged.
!  (in)      GOLD_dom - The GOLD_domain_type containing the mpp_domain needed to
!                      determine where data should be sent.
!  (in)      sideflag - An optional integer indicating which directions the
!                       data should be sent.  It is TO_ALL or the sum of any of
!                       TO_EAST, TO_WEST, TO_NORTH, and TO_SOUTH.  For example,
!                       TO_EAST sends the data to the processor to the east, so
!                       the halos on the western side are filled.  TO_ALL is
!                       the default if sideflag is omitted.
!  (in)      position - An optional argument indicating the position.  This is
!                       may be CORNER, but is CENTER by default.
!  (in)      complete - An optional argument indicating whether the halo updates
!                       should be initiated immediately or wait for second 
!                       pass_..._start call.  Omitting complete is the same as 
!                       setting complete to .true.
!  (return value) - The integer index for this update.
  integer :: dirflag

  dirflag = To_All ! 60
  if (PRESENT(sideflag)) then ; if (sideflag > 0) dirflag = sideflag ; endif

#ifdef FORBID_NONBLOCKING_UPDATES
  call GOLD_error(FATAL, "pass_var_start: This call is not yet available.")
  pass_var_start_3d = -1
#else
  pass_var_start_3d = mpp_start_update_domains(array, GOLD_dom%mpp_domain, &
                          flags=dirflag, position=position)
#endif
end function pass_var_start_3d

subroutine pass_var_complete_2d(id_update, array, GOLD_dom, sideflag, position)
  integer,                intent(in)    :: id_update
  real, dimension(:,:),   intent(inout) :: array
  type(GOLD_domain_type), intent(inout) :: GOLD_dom
  integer,      optional, intent(in)    :: sideflag
  integer,      optional, intent(in)    :: position
  integer :: pass_var_start_2d
! Arguments: id_update - The integer id of this update which has been returned
!                        from a previous call to pass_var_start.
!  (inout)   array - The array which is having its halos points exchanged.
!  (in)      GOLD_dom - The GOLD_domain_type containing the mpp_domain needed to
!                      determine where data should be sent.
!  (in)      sideflag - An optional integer indicating which directions the
!                       data should be sent.  It is TO_ALL or the sum of any of
!                       TO_EAST, TO_WEST, TO_NORTH, and TO_SOUTH.  For example,
!                       TO_EAST sends the data to the processor to the east, so
!                       the halos on the western side are filled.  TO_ALL is
!                       the default if sideflag is omitted.
!  (in)      position - An optional argument indicating the position.  This is
!                       may be CORNER, but is CENTER by default.
  integer :: dirflag

  dirflag = To_All ! 60
  if (PRESENT(sideflag)) then ; if (sideflag > 0) dirflag = sideflag ; endif

#ifdef FORBID_NONBLOCKING_UPDATES
  call GOLD_error(FATAL, "pass_var_complete: This call is not yet available.")
#else
  call mpp_complete_update_domains(id_update, array, GOLD_dom%mpp_domain, &
                                   flags=dirflag, position=position)
#endif
end subroutine pass_var_complete_2d

subroutine pass_var_complete_3d(id_update, array, GOLD_dom, sideflag, position)
  integer,                intent(in)    :: id_update
  real, dimension(:,:,:), intent(inout) :: array
  type(GOLD_domain_type), intent(inout) :: GOLD_dom
  integer,      optional, intent(in)    :: sideflag
  integer,      optional, intent(in)    :: position
! Arguments: id_update - The integer id of this update which has been returned
!                        from a previous call to pass_var_start.
!  (inout)   array - The array which is having its halos points exchanged.
!  (in)      GOLD_dom - The GOLD_domain_type containing the mpp_domain needed to
!                      determine where data should be sent.
!  (in)      sideflag - An optional integer indicating which directions the
!                       data should be sent.  It is TO_ALL or the sum of any of
!                       TO_EAST, TO_WEST, TO_NORTH, and TO_SOUTH.  For example,
!                       TO_EAST sends the data to the processor to the east, so
!                       the halos on the western side are filled.  TO_ALL is
!                       the default if sideflag is omitted.
!  (in)      position - An optional argument indicating the position.  This is
!                       may be CORNER, but is CENTER by default.
  integer :: dirflag

  dirflag = To_All ! 60
  if (PRESENT(sideflag)) then ; if (sideflag > 0) dirflag = sideflag ; endif

#ifdef FORBID_NONBLOCKING_UPDATES
  call GOLD_error(FATAL, "pass_var_complete: This call is not yet available.")
#else
  call mpp_complete_update_domains(id_update, array, GOLD_dom%mpp_domain, &
                                   flags=dirflag, position=position)
#endif
end subroutine pass_var_complete_3d


subroutine pass_vector_2d(u_cmpt, v_cmpt, GOLD_dom, direction, stagger, complete)
  real, dimension(:,:),  intent(inout) :: u_cmpt, v_cmpt
  type(GOLD_domain_type), intent(inout) :: GOLD_dom
  integer,     optional, intent(in)    :: direction
  integer,     optional, intent(in)    :: stagger
  logical,     optional, intent(in)    :: complete
! Arguments: u_cmpt - The nominal zonal (u) component of the vector pair which
!                     is having its halos points exchanged.
!  (inout)   v_cmpt - The nominal meridional (v) component of the vector pair
!                     which is having its halos points exchanged. 
!  (in)      GOLD_dom - The GOLD_domain_type containing the mpp_domain needed to
!                      determine where data should be sent.
!  (in)      direction - An optional integer indicating which directions the
!                        data should be sent.  It is TO_ALL or the sum of any of
!                        TO_EAST, TO_WEST, TO_NORTH, and TO_SOUTH, possibly
!                        plus SCALAR_PAIR if these are paired non-directional
!                        scalars discretized at the typical vector component
!                        locations.  For example, TO_EAST sends the data to the
!                        processor to the east, so the halos on the western
!                        side are filled.  TO_ALL is the default if omitted.
!  (in)      stagger - An optional flag, which may be one of A_GRID, BGRID_NE,
!                      or CGRID_NE, indicating where the two components of the
!                      vector are discretized.  Omitting stagger is the same as
!                      setting it to CGRID_NE.
!  (in)      complete - An optional argument indicating whether the halo updates
!                       should be completed before progress resumes.  Omitting
!                       complete is the same as setting complete to .true.

  integer :: stagger_local
  integer :: dirflag
  logical :: block_til_complete

  stagger_local = CGRID_NE ! Default value for type of grid
  if (present(stagger)) stagger_local = stagger

  dirflag = To_All ! 60
  if (PRESENT(direction)) then ; if (direction > 0) dirflag = direction ; endif
  block_til_complete = .true.
  if (present(complete)) block_til_complete = complete

  call mpp_update_domains(u_cmpt, v_cmpt, GOLD_dom%mpp_domain, flags=dirflag, &
                          gridtype=stagger_local, complete = block_til_complete)
end subroutine pass_vector_2d

subroutine pass_vector_3d(u_cmpt, v_cmpt, GOLD_dom, direction, stagger, complete)
  real, dimension(:,:,:), intent(inout) :: u_cmpt, v_cmpt
  type(GOLD_domain_type),  intent(inout) :: GOLD_dom
  integer,      optional, intent(in)    :: direction
  integer,      optional, intent(in)    :: stagger
  logical,      optional, intent(in)    :: complete
! Arguments: u_cmpt - The nominal zonal (u) component of the vector pair which
!                     is having its halos points exchanged.
!  (inout)   v_cmpt - The nominal meridional (v) component of the vector pair
!                     which is having its halos points exchanged. 
!  (in)      GOLD_dom - The GOLD_domain_type containing the mpp_domain needed to
!                      determine where data should be sent.
!  (in)      direction - An optional integer indicating which directions the
!                        data should be sent.  It is TO_ALL or the sum of any of
!                        TO_EAST, TO_WEST, TO_NORTH, and TO_SOUTH, possibly
!                        plus SCALAR_PAIR if these are paired non-directional
!                        scalars discretized at the typical vector component
!                        locations.  For example, TO_EAST sends the data to the
!                        processor to the east, so the halos on the western
!                        side are filled.  TO_ALL is the default if omitted.
!  (in)      stagger - An optional flag, which may be one of A_GRID, BGRID_NE,
!                      or CGRID_NE, indicating where the two components of the
!                      vector are discretized.  Omitting stagger is the same as
!                      setting it to CGRID_NE.
!  (in)      complete - An optional argument indicating whether the halo updates
!                       should be completed before progress resumes.  Omitting
!                       complete is the same as setting complete to .true.

  integer :: stagger_local
  integer :: dirflag
  logical :: block_til_complete

  stagger_local = CGRID_NE ! Default value for type of grid
  if (present(stagger)) stagger_local = stagger

  dirflag = To_All ! 60
  if (PRESENT(direction)) then ; if (direction > 0) dirflag = direction ; endif
  block_til_complete = .true.
  if (present(complete)) block_til_complete = complete

  call mpp_update_domains(u_cmpt, v_cmpt, GOLD_dom%mpp_domain, flags=dirflag, &
                          gridtype=stagger_local, complete = block_til_complete)
end subroutine pass_vector_3d

function pass_vector_start_2d(u_cmpt, v_cmpt, GOLD_dom, direction, stagger, complete)
  real, dimension(:,:),   intent(inout) :: u_cmpt, v_cmpt
  type(GOLD_domain_type), intent(inout) :: GOLD_dom
  integer,      optional, intent(in)    :: direction
  integer,      optional, intent(in)    :: stagger
  logical,      optional, intent(in)    :: complete
  integer                               :: pass_vector_start_2d
! Arguments: u_cmpt - The nominal zonal (u) component of the vector pair which
!                     is having its halos points exchanged.
!  (inout)   v_cmpt - The nominal meridional (v) component of the vector pair
!                     which is having its halos points exchanged. 
!  (in)      GOLD_dom - The GOLD_domain_type containing the mpp_domain needed to
!                      determine where data should be sent.
!  (in)      direction - An optional integer indicating which directions the
!                        data should be sent.  It is TO_ALL or the sum of any of
!                        TO_EAST, TO_WEST, TO_NORTH, and TO_SOUTH, possibly
!                        plus SCALAR_PAIR if these are paired non-directional
!                        scalars discretized at the typical vector component
!                        locations.  For example, TO_EAST sends the data to the
!                        processor to the east, so the halos on the western
!                        side are filled.  TO_ALL is the default if omitted.
!  (in)      stagger - An optional flag, which may be one of A_GRID, BGRID_NE,
!                      or CGRID_NE, indicating where the two components of the
!                      vector are discretized.  Omitting stagger is the same as
!                      setting it to CGRID_NE.
!  (in)      complete - An optional argument indicating whether the halo updates
!                       should be initiated immediately or wait for second 
!                       pass_..._start call.  Omitting complete is the same as 
!                       setting complete to .true.
!  (return value) - The integer index for this update.
  integer :: stagger_local
  integer :: dirflag

  stagger_local = CGRID_NE ! Default value for type of grid
  if (present(stagger)) stagger_local = stagger

  dirflag = To_All ! 60
  if (PRESENT(direction)) then ; if (direction > 0) dirflag = direction ; endif

#ifdef FORBID_NONBLOCKING_UPDATES
  call GOLD_error(FATAL, "pass_vector_start: This call is not yet available.")
  pass_vector_start_2d = -1
#else
  pass_vector_start_2d = mpp_start_update_domains(u_cmpt, v_cmpt, &
      GOLD_dom%mpp_domain, flags=dirflag, gridtype=stagger_local)
#endif

end function pass_vector_start_2d

function pass_vector_start_3d(u_cmpt, v_cmpt, GOLD_dom, direction, stagger, complete)
  real, dimension(:,:,:), intent(inout) :: u_cmpt, v_cmpt
  type(GOLD_domain_type), intent(inout) :: GOLD_dom
  integer,      optional, intent(in)    :: direction
  integer,      optional, intent(in)    :: stagger
  logical,      optional, intent(in)    :: complete
  integer                               :: pass_vector_start_3d
! Arguments: u_cmpt - The nominal zonal (u) component of the vector pair which
!                     is having its halos points exchanged.
!  (inout)   v_cmpt - The nominal meridional (v) component of the vector pair
!                     which is having its halos points exchanged. 
!  (in)      GOLD_dom - The GOLD_domain_type containing the mpp_domain needed to
!                      determine where data should be sent.
!  (in)      direction - An optional integer indicating which directions the
!                        data should be sent.  It is TO_ALL or the sum of any of
!                        TO_EAST, TO_WEST, TO_NORTH, and TO_SOUTH, possibly
!                        plus SCALAR_PAIR if these are paired non-directional
!                        scalars discretized at the typical vector component
!                        locations.  For example, TO_EAST sends the data to the
!                        processor to the east, so the halos on the western
!                        side are filled.  TO_ALL is the default if omitted.
!  (in)      stagger - An optional flag, which may be one of A_GRID, BGRID_NE,
!                      or CGRID_NE, indicating where the two components of the
!                      vector are discretized.  Omitting stagger is the same as
!                      setting it to CGRID_NE.
!  (in)      complete - An optional argument indicating whether the halo updates
!                       should be initiated immediately or wait for second 
!                       pass_..._start call.  Omitting complete is the same as 
!                       setting complete to .true.
!  (return value) - The integer index for this update.
  integer :: stagger_local
  integer :: dirflag

  stagger_local = CGRID_NE ! Default value for type of grid
  if (present(stagger)) stagger_local = stagger

  dirflag = To_All ! 60
  if (PRESENT(direction)) then ; if (direction > 0) dirflag = direction ; endif

#ifdef FORBID_NONBLOCKING_UPDATES
  call GOLD_error(FATAL, "pass_vector_start: This call is not yet available.")
  pass_vector_start_3d = -1
#else
  pass_vector_start_3d = mpp_start_update_domains(u_cmpt, v_cmpt, &
      GOLD_dom%mpp_domain, flags=dirflag, gridtype=stagger_local)
#endif

end function pass_vector_start_3d

subroutine pass_vector_complete_2d(id_update, u_cmpt, v_cmpt, GOLD_dom, direction, stagger)
  integer,                intent(in)    :: id_update
  real, dimension(:,:),   intent(inout) :: u_cmpt, v_cmpt
  type(GOLD_domain_type), intent(inout) :: GOLD_dom
  integer,      optional, intent(in)    :: direction
  integer,      optional, intent(in)    :: stagger
  integer                               :: pass_vector_start_2d
! Arguments: id_update - The integer id of this update which has been returned
!                        from a previous call to pass_var_start.
!  (inout)   u_cmpt - The nominal zonal (u) component of the vector pair which
!                     is having its halos points exchanged.
!  (inout)   v_cmpt - The nominal meridional (v) component of the vector pair
!                     which is having its halos points exchanged. 
!  (in)      GOLD_dom - The GOLD_domain_type containing the mpp_domain needed to
!                      determine where data should be sent.
!  (in)      direction - An optional integer indicating which directions the
!                        data should be sent.  It is TO_ALL or the sum of any of
!                        TO_EAST, TO_WEST, TO_NORTH, and TO_SOUTH, possibly
!                        plus SCALAR_PAIR if these are paired non-directional
!                        scalars discretized at the typical vector component
!                        locations.  For example, TO_EAST sends the data to the
!                        processor to the east, so the halos on the western
!                        side are filled.  TO_ALL is the default if omitted.
!  (in)      stagger - An optional flag, which may be one of A_GRID, BGRID_NE,
!                      or CGRID_NE, indicating where the two components of the
!                      vector are discretized.  Omitting stagger is the same as
!                      setting it to CGRID_NE.
!  (return value) - The integer index for this update.
  integer :: stagger_local
  integer :: dirflag

  stagger_local = CGRID_NE ! Default value for type of grid
  if (present(stagger)) stagger_local = stagger

  dirflag = To_All ! 60
  if (PRESENT(direction)) then ; if (direction > 0) dirflag = direction ; endif

#ifdef FORBID_NONBLOCKING_UPDATES
  call GOLD_error(FATAL, "pass_vector_complete: This call is not yet available.")
#else
  call mpp_complete_update_domains(id_update, u_cmpt, v_cmpt, &
           GOLD_dom%mpp_domain, flags=dirflag, gridtype=stagger_local)
#endif

end subroutine pass_vector_complete_2d

subroutine pass_vector_complete_3d(id_update, u_cmpt, v_cmpt, GOLD_dom, direction, stagger)
  integer,                intent(in)    :: id_update
  real, dimension(:,:,:), intent(inout) :: u_cmpt, v_cmpt
  type(GOLD_domain_type), intent(inout) :: GOLD_dom
  integer,      optional, intent(in)    :: direction
  integer,      optional, intent(in)    :: stagger
  integer                               :: pass_vector_start_2d
! Arguments: id_update - The integer id of this update which has been returned
!                        from a previous call to pass_var_start.
!  (inout)   u_cmpt - The nominal zonal (u) component of the vector pair which
!                     is having its halos points exchanged.
!  (inout)   v_cmpt - The nominal meridional (v) component of the vector pair
!                     which is having its halos points exchanged. 
!  (in)      GOLD_dom - The GOLD_domain_type containing the mpp_domain needed to
!                      determine where data should be sent.
!  (in)      direction - An optional integer indicating which directions the
!                        data should be sent.  It is TO_ALL or the sum of any of
!                        TO_EAST, TO_WEST, TO_NORTH, and TO_SOUTH, possibly
!                        plus SCALAR_PAIR if these are paired non-directional
!                        scalars discretized at the typical vector component
!                        locations.  For example, TO_EAST sends the data to the
!                        processor to the east, so the halos on the western
!                        side are filled.  TO_ALL is the default if omitted.
!  (in)      stagger - An optional flag, which may be one of A_GRID, BGRID_NE,
!                      or CGRID_NE, indicating where the two components of the
!                      vector are discretized.  Omitting stagger is the same as
!                      setting it to CGRID_NE.
!  (return value) - The integer index for this update.
  integer :: stagger_local
  integer :: dirflag

  stagger_local = CGRID_NE ! Default value for type of grid
  if (present(stagger)) stagger_local = stagger

  dirflag = To_All ! 60
  if (PRESENT(direction)) then ; if (direction > 0) dirflag = direction ; endif

#ifdef FORBID_NONBLOCKING_UPDATES
  call GOLD_error(FATAL, "pass_vector_complete: This call is not yet available.")
#else
  call mpp_complete_update_domains(id_update, u_cmpt, v_cmpt, &
           GOLD_dom%mpp_domain, flags=dirflag, gridtype=stagger_local)
#endif

end subroutine pass_vector_complete_3d

subroutine GOLD_domains_init(GOLD_dom, param_file, min_halo, symmetric)
  type(GOLD_domain_type),          pointer       :: GOLD_dom
  type(param_file_type),           intent(in)    :: param_file
  integer, dimension(2), optional, intent(inout) :: min_halo
  logical, optional,               intent(in)    :: symmetric
! Arguments: GOLD_dom - A pointer to the GOLD_domain_type being defined here.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in,opt)  min_halo - If present, this sets the minimum halo size for this
!                       domain in the x- and y- directions, and returns the
!                       actual halo size used.
!  (in,opt)  symmetric - If present, this specified whether this domain is
!                        symmetric, regardless of whether SYMMETRIC_MEMORY is
!                        defined.

  integer, dimension(2) :: layout = (/ 1, 1 /)
  integer, dimension(2) :: io_layout = (/ 0, 0 /)
  integer, dimension(4) :: global_indices
  integer :: x_halo, y_halo
  integer :: pe, proc_used
  integer :: isc,iec,jsc,jec ! The bounding indices of the computational domain.
  integer :: X_FLAGS, Y_FLAGS
  integer :: i, xsiz, ysiz
  logical :: reentrant_x, reentrant_y, tripolar_N, is_static
  character(len=200) :: mesg
  character(len=8) :: char_xsiz, char_ysiz, char_NXTOT, char_NYTOT
  character(len=128) :: version = '$Id$'
  character(len=128) :: tagname = '$Name$'
  character(len=40)  :: mod ! This module's name.

  if (.not.associated(GOLD_dom)) then
    allocate(GOLD_dom)
    allocate(GOLD_dom%mpp_domain)
  endif

  pe = PE_here()
  proc_used = num_PEs()

  if (present(symmetric)) then
    GOLD_dom%symmetric = symmetric
    if (symmetric) then ; mod = "GOLD_domains symmetric"
    else ; mod = "GOLD_domains non-symmetric" ; endif
  else
    mod = "GOLD_domains"
#ifdef SYMMETRIC_MEMORY
    GOLD_dom%symmetric = .true.
#else
    GOLD_dom%symmetric = .false.
#endif
  endif
  if (present(min_halo)) mod = trim(mod)//" min_halo"
  
  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mod, version, tagname)
  call get_param(param_file, mod, "REENTRANT_X", reentrant_x, &
                 "If true, the domain is zonally reentrant.", default=.true.)
  call get_param(param_file, mod, "REENTRANT_Y", reentrant_y, &
                 "If true, the domain is meridionally reentrant.", &
                 default=.false.)
  call get_param(param_file, mod, "TRIPOLAR_N", tripolar_N, &
                 "Use tripolar connectivity at the northern edge of the \n"//&
                 "domain.  With TRIPOLAR_N, NXTOT must be even.", &
                 default=.false.)
  
  call log_param(param_file, mod, "SYMMETRIC_MEMORY", GOLD_dom%symmetric, &
                 "If defined, the velocity point data domain includes \n"//&
                 "every face of the thickness points. In other words, \n"//&
                 "some arrays are larger than others, depending on where \n"//&
                 "they are on the staggered grid.  Also, the starting \n"//&
                 "index of the velocity-point arrays is usually 0, not 1. \n"//&
                 "This can only be set at compile time.")
  call get_param(param_file, mod, "NONBLOCKING_UPDATES", GOLD_dom%nonblocking_updates, &
                 "If true, non-blocking halo updates may be used.", &
                 default=.false.)
#ifdef FORBID_NONBLOCKING_UPDATES
  if (GOLD_dom%nonblocking_updates) call GOLD_error(FATAL, "GOLD_domains_init: "//&
         "Non-blocking halo updates are not yet permitted.  Use #undef NONBLOCKING_UPDATES.")
#endif
  is_static = .false.
#ifdef STATIC_MEMORY
  is_static = .true.
#endif
  call log_param(param_file, mod, "STATIC_MEMORY", is_static, &
                 "If STATIC_MEMORY is defined, the principle variables \n"//&
                 "will have sizes that are statically determined at \n"//&
                 "compile time.  Otherwise the sizes are not determined \n"//&
                 "until run time. The STATIC option is substantially \n"//&
                 "faster, but does not allow the PE count to be changed \n"//&
                 "at run time.  This can only be set at compile time.")

  call get_param(param_file, mod, "NXHALO", GOLD_dom%nx_halo, &
                 "The number of halo points on each side in the \n"//&
                 "x-direction.  With STATIC_MEMORY this is set as NX_HALO \n"//&
                 "in GOLD_memory.h at compile time; without STATIC_MEMORY \n"//&
                 "the default is NX_HALO in GOLD_memory.h.", default=NX_HALO)
  call get_param(param_file, mod, "NYHALO", GOLD_dom%ny_halo, &
                 "The number of halo points on each side in the \n"//&
                 "y-direction.  With STATIC_MEMORY this is set as NY_HALO \n"//&
                 "in GOLD_memory.h at compile time; without STATIC_MEMORY \n"//&
                 "the default is NY_HALO in GOLD_memory.h.", default=NY_HALO)
  if (present(min_halo)) then
    GOLD_dom%nx_halo = max(GOLD_dom%nx_halo, min_halo(1))
    min_halo(1) = GOLD_dom%nx_halo
    GOLD_dom%ny_halo = max(GOLD_dom%ny_halo, min_halo(2))
    min_halo(2) = GOLD_dom%ny_halo
    call log_param(param_file, mod, "NXHALO min_halo", GOLD_dom%nx_halo)
    call log_param(param_file, mod, "NYHALO min_halo", GOLD_dom%nx_halo)
  endif
#ifdef STATIC_MEMORY
  call get_param(param_file, mod, "NXTOT", GOLD_dom%nxtot, &
                 "The total number of thickness grid points in the \n"//&
                 "x-direction in the physical domain. With STATIC_MEMORY \n"//&
                 "this is set in GOLD_memory.h at compile time.", default=NXTOT)
  call get_param(param_file, mod, "NYTOT", GOLD_dom%nytot, &
                 "The total number of thickness grid points in the \n"//&
                 "x-direction in the physical domain. With STATIC_MEMORY \n"//&
                 "this is set in GOLD_memory.h at compile time.", default=NYTOT)
  if (GOLD_dom%nxtot /= NXTOT) call GOLD_error(FATAL,"GOLD_domains_init: " // &
   "static mismatch for NXTOT domain size. Header file does not match input namelist")
  if (GOLD_dom%nytot /= NYTOT) call GOLD_error(FATAL,"GOLD_domains_init: " // &
   "static mismatch for NYTOT domain size. Header file does not match input namelist")

  if (.not.present(min_halo)) then
    if (GOLD_dom%nx_halo /= NX_HALO) call GOLD_error(FATAL,"GOLD_domains_init: " // &
           "static mismatch for NXHALO domain size")
    if (GOLD_dom%ny_halo /= NY_HALO) call GOLD_error(FATAL,"GOLD_domains_init: " // &
           "static mismatch for NYHALO domain size")
  endif
#else
  call get_param(param_file, mod, "NXTOT", GOLD_dom%nxtot, &
                 "The total number of thickness grid points in the \n"//&
                 "x-direction in the physical domain. With STATIC_MEMORY \n"//&
                 "this is set in GOLD_memory.h at compile time.", &
                 fail_if_missing=.true.)
  call get_param(param_file, mod, "NYTOT", GOLD_dom%nytot, &
                 "The total number of thickness grid points in the \n"//&
                 "x-direction in the physical domain. With STATIC_MEMORY \n"//&
                 "this is set in GOLD_memory.h at compile time.", &
                 fail_if_missing=.true.)
#endif
  x_halo = GOLD_dom%nx_halo
  y_halo = GOLD_dom%ny_halo

  global_indices(1) = 1
  global_indices(2) = GOLD_dom%nxtot
  global_indices(3) = 1
  global_indices(4) = GOLD_dom%nytot

#ifdef STATIC_MEMORY
  layout(1) = NXPROC ; layout(2) = NYPROC
#else
  call mpp_define_layout(global_indices, proc_used, layout)
  call read_param(param_file,"NXPROC",layout(1))
  call read_param(param_file,"NYPROC",layout(2))
  if (layout(1)*layout(2) /= proc_used) then
    write(mesg,'("GOLD_domains_init: The product of the two components of layout, ", &
          &      2i4,", is not the number of PEs used, ",i5,".")') &
          layout(1),layout(2),proc_used
    call GOLD_error(FATAL, mesg)
  endif
#endif
  call log_param(param_file, mod, "NXPROC", layout(1), &
                 "The number of processors in the x-direction. With \n"//&
                 "STATIC_MEMORY this is set in GOLD_memory.h at compile time.")
  call log_param(param_file, mod, "NYPROC", layout(2), &
                 "The number of processors in the x-direction. With \n"//&
                 "STATIC_MEMORY this is set in GOLD_memory.h at compile time.")
!  write(*,*) 'layout is now ',layout, global_indices

!if (layout(1) < NXPROC) &
!  call GOLD_error(FATAL,"GOLD_domains: layout in zonal axis has too few processors.")
!if (layout(2) < NYPROC) &
!  call GOLD_error(FATAL,"GOLD_domains: layout in meridional axis has too few processors.")

  !   Set up the I/O lay-out, and check that it uses an even multiple of the
  ! number of PEs in each direction.
  io_layout(:) = (/ 0, 0 /)
  call get_param(param_file, mod, "NXPROC_IO", io_layout(1), &
                 "The number of processors used for I/O in the \n"//&
                 "x-direction, or 0 to equal NXPROC.", default=0)
  call get_param(param_file, mod, "NYPROC_IO", io_layout(2), &
                 "The number of processors used for I/O in the \n"//&
                 "y-direction, or 0 to equal NYPROC.", default=0)
  if (io_layout(1) < 0) then
    write(mesg,'("GOLD_domains_init: NXPROC_IO = ",i4,".  Negative values of "//&
         &" of NXPROC_IO are not allowed.")') io_layout(1)
    call GOLD_error(FATAL, mesg)
  elseif (io_layout(1) > 0) then ; if (modulo(layout(1), io_layout(1)) /= 0) then
    write(mesg,'("GOLD_domains_init: The x-direction I/O-layout, NXPROC_IO=",i4, &
         &", does not evenly divide the x-direction layout, NXPROC=,",i4,".")') &
          io_layout(1),layout(1)
    call GOLD_error(FATAL, mesg)
  endif ; endif
  
  if (io_layout(2) < 0) then
    write(mesg,'("GOLD_domains_init: NYPROC_IO = ",i4,".  Negative values of "//&
         &" of NYPROC_IO are not allowed.")') io_layout(2)
    call GOLD_error(FATAL, mesg)
  elseif (io_layout(2) /= 0) then ; if (modulo(layout(2), io_layout(2)) /= 0) then
    write(mesg,'("GOLD_domains_init: The y-direction I/O-layout, NYPROC_IO=",i4, &
         &", does not evenly divide the y-direction layout, NYPROC=,",i4,".")') &
          io_layout(2),layout(2)
    call GOLD_error(FATAL, mesg)
  endif ; endif

  X_FLAGS = 0 ; Y_FLAGS = 0
  if (reentrant_x) X_FLAGS = CYCLIC_GLOBAL_DOMAIN
  if (reentrant_y) Y_FLAGS = CYCLIC_GLOBAL_DOMAIN
  if (tripolar_N) then
    Y_FLAGS = FOLD_NORTH_EDGE
    if (reentrant_y) call GOLD_error(FATAL,"GOLD_domains: "// &
      "TRIPOLAR_N and REENTRANT_Y may not be defined together.")
  endif
  
  call GOLD_define_domain((/1+x_halo,GOLD_dom%nxtot+x_halo,1+y_halo, &
                 GOLD_dom%nytot+y_halo/), layout, GOLD_dom%mpp_domain, &
                 xflags=X_FLAGS, yflags=Y_FLAGS, xhalo=x_halo, yhalo=y_halo, &
                 symmetry = GOLD_dom%symmetric, name="GOLD")

  if ((io_layout(1) + io_layout(2) > 0)) then
    call GOLD_define_io_domain(GOLD_dom%mpp_domain, io_layout)
  endif

! Save the extra data for creating other domains of different resolution that overlay this domain
  GOLD_dom%X_FLAGS = X_FLAGS
  GOLD_dom%Y_FLAGS = Y_FLAGS
  GOLD_dom%layout = layout
  GOLD_dom%io_layout = io_layout
  GOLD_dom%use_io_layout = (io_layout(1) + io_layout(2) > 0)

#ifdef STATIC_MEMORY
!   A requirement of equal sized compute domains is necessary when STATIC_MEMORY
! is used.
  call mpp_get_compute_domain(GOLD_dom%mpp_domain,isc,iec,jsc,jec)
  xsiz = iec - isc + 1
  ysiz = jec - jsc + 1
  if (xsiz*NXPROC /= GOLD_dom%nxtot .OR. ysiz*NYPROC /= GOLD_dom%nytot) then
     write( char_xsiz,'(i4)' ) NXPROC
     write( char_ysiz,'(i4)' ) NYPROC
     write( char_NXTOT,'(i4)' ) GOLD_dom%nxtot
     write( char_NYTOT,'(i4)' ) GOLD_dom%nytot
     call GOLD_error(WARNING,'GOLD_domains: Processor decomposition (NXPROC,NYPROC) = (' &
         //trim(char_xsiz)//','//trim(char_ysiz)// &
         ') does not evenly divide size set by preprocessor macro ('&
         //trim(char_NXTOT)//','//trim(char_NYTOT)// '). ')
     call GOLD_error(FATAL,'GOLD_domains:  #undef STATIC_MEMORY in GOLD_memory.h to use &
         &dynamic allocation, or change processor decomposition to evenly divide the domain.')
  endif
#endif

end subroutine GOLD_domains_init

subroutine get_domain_extent(Domain, isc, iec, jsc, jec, isd, ied, jsd, jed, &
                             isg, ieg, jsg, jeg, idg_offset, jdg_offset, symmetric)
  type(GOLD_domain_type), intent(in) :: Domain
  integer, intent(out) :: isc, iec, jsc, jec
  integer, intent(out) :: isd, ied, jsd, jed
  integer, intent(out) :: isg, ieg, jsg, jeg
  integer, intent(out) :: idg_offset, jdg_offset
  logical, intent(out) :: symmetric
! Arguments: Domain - The GOLD_domain_type from which the indices are extracted.
!  (out)     isc, iec, jsc, jec - the start & end indices of the
!                                 computational domain.
!  (out)     isd, ied, jsd, jed - the start & end indices of the data domain.
!  (out)     isg, ieg, jsg, jeg - the start & end indices of the global domain.
!  (out)     idg_offset, jdg_offset - the offset between the corresponding
!                                     global and data index spaces.
!  (out)     symmetric - true if symmetric memory is used.

  call mpp_get_compute_domain(Domain%mpp_domain, isc, iec, jsc, jec)
  call mpp_get_data_domain(Domain%mpp_domain, isd, ied, jsd, jed)
  call mpp_get_global_domain(Domain%mpp_domain, isg, ieg, jsg, jeg)

  ! This code institutes the GOLD convention that local array indices start at 1.
  idg_offset = isd-1 ; jdg_offset = jsd-1
  isc = isc-isd+1 ; iec = iec-isd+1 ; jsc = jsc-jsd+1 ; jec = jec-jsd+1
  ied = ied-isd+1 ; jed = jed-jsd+1
  isd = 1 ; jsd = 1
  symmetric = Domain%symmetric

end subroutine get_domain_extent

end module GOLD_domains
